Dense-dilute factorization for a class of stochastic processes and for high energy QCD 
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Stochastic processes described by evolution equations in the universality class of the FKPP equa- 
tion may be approximately factorized into a linear stochastic part and a nonlinear deterministic 
part. We prove this factorization on a model with no spatial dimensions and we illustrate it nu- 
merically on a one-dimensional toy model that possesses some of the main features of high energy 
QCD evolution. We explain how this procedure may be applied to QCD amplitudes, by combining 
Salam's Monte-Carlo implementation of the dipole model and a numerical solution of the Balitsky- 
Kovchegov equation. 
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I. INTRODUCTION 

High energy scattering in QCD was recently shown 
to be essentially similar to a reaction-diffusion process 
[HQ- To understand this correspondence, one needs to 
introduce an (unphysical) forward elastic scattering am- 
plitude T(Y,r), that encodes the interaction probability 
for given parton or field configurations of the incoming 
hadrons, at rapidity Y, and for a characteristic trans- 
verse distance scale r. The physical scattering amplitude 
A(Y, r) is the average of T(Y, r) over all possible realiza- 
tions of the fields or partons. Although T is not measur- 
able experimentally, it is a quantity that has also to be 
introduced when one wants to write a Monte-Carlo event 
generator. 

The equation that governs the rapidity evolution of 
T at fixed impact-parameter belongs to the universal- 
ity class of the stochastic Fishcr-Kolmogorov-Petrovsky- 
Piscounov (FKPP) equation Q, or equivalently, of 
the Reggeon field theory equation. It can be seen 
as a stochastic extension of the well-known Balitsky- 
Kovchegov (BK) equation [J, [j| , that describes a peculiar 
limit of high energy QCD. It is convenient 1 to go to mo- 
mentum space using the transformation [f| 



f(Y, k) = 



Akr 



27T7- 2 



T(Y,r) 



(1) 



The evolution of T can be written, for example, in the 
form 



d &Y f = x(-d lnk2 )f-f 2 



2Tv, 



(2) 



where a = a s N c /ir, x is a representation of the BFKL 
kernel @ in momentum space and v is a noise of zero 
mean that varies randomly by typically one unit when 
aY or lnfc 2 are changed by one unit. A is the average of 
T with respect to the noise v. 

This new approach to high energy scattering in QCD, 
that makes use of ideas and tools of statistical physics, 



The reason for introducing this transformation is that the non- 
linearity present in the Balitsky-Kovchegov equation simplifies 
greatly in momentum space. 



has already inspired a number of theoretical works. 
Eq. itself was also subsequently discussed in Ref. Q 
and studied numerically (with a Gaussian noise v) in 
Ref. [|[. Many developments from different perspectives 
have followed: In particular, the connection to Reggeon 
field theory was investi gate d Q, a QCD derivation of 
Eq. ([2]) was searched for [10| , and different formalisms to 
describe the same physics were proposed fllj [. 

Eq. ^ may be interpreted in the following way: It 
describes the evolution in time t — aY of a fraction T of 
N = 1/ev 2 . particles (gluons) per unit of x = Infc 2 that 
multiply and diffuse in x through the branching diffusion 
kernel \, and that recombine through the nonlinear term. 
Eq. ([2]) is however not an exact equation of QCD, and 
should by no way be interpreted as such. It is rather a 
synthetic form of writing two exactly known limits: the 
large parton density limit, in which the evolution of T is 
given by the Balitsky-Kovchegov equation 



d aY T = x (-d lnk 2)T-T 2 , 



(3) 



and the low density limit, 
stochastic equation 



d &Y T = x(~d lnk 2)T + a s V2Tv. 



represented by the linear 



(4) 



Let us briefly recall a bit more precisely the QCD con- 
tent of Eq. @ (We refer the reader to the original pa- 
pers [H, Q for details). In the low density region of phase 
space (equivalently, in the region in which the amplitude 
is small, T,T <C 1), it is useful to view rapidity evo- 
lution in the framework of the color dipole model [l2 |. 
One goes to the reference frame of one of the interacting 
hadrons, that we assume to be a quark- antiquark dipole 
for simplicity and that we will call the probe. The sec- 
ond hadron carries all the rapidity and QCD evolution 
and will be referred to as the target. In the large- N c 
limit, the Fock state of the target may be represented by 
a set of color dipoles. Each of these dipoles, character- 
ized by a two-dimensional position vector xqi — xq — x\ 
in transverse coordinate space, is spanned by two gluons 
respectively sitting at positions xo and x\. The model 
is defined by the stochastic branching of the dipoles that 
occurs as the target is boosted to larger rapidities: In a 
step dY of rapidity, each dipole present at rapidity Y has 



2 



the probability 



adY 



d 2 x 2 



'01 



2n 



X 02 X 12 



(5) 



to be replaced by two new dipoles of respective sizes X02 
and x\2- 

When nonlinear effects can be neglected, the relation- 
ship between the number of dipoles n(Y, Xqi) and the am- 
plitude T(Y, r i) for the scattering of the probe dipole of 
size roi off a random configuration of the target reads, in 
appropriate normalizations (see e.g. [Hi), 



d 2 x d 2 xi 2 |r - a;i| 2 |ri 
x / — — In 



2tt 2tt 



|r - x \ 2 \ri - xi\ 



■n(Y,x 01 ). (6) 



This relationship turns out to be approximately local in 
dipole sizes and impact parameter: T(Y, roi) is of order 
a 2 times the number of dipoles in the considered config- 
uration of the target whose sizes are in a bin of width 1 
(on a logarithmic scale) centered on |roi|, and which sit 
within a distance |tqi | of the probe in impact parameter. 

By converting the splitting process ((5|) into an evolu- 
tion equation for T with the help of Eq. (J6j) , by transform- 
ing to momentum space using Eq. (fT]), and by averaging 
over the angle in the transverse plane, one would get the 
linear terms in Eq. together with an appropriate noise 
v that would encode the fluctuations in the dipole num- 
ber induced by the stochastic branching process ([5]) : this 
is precisely Eq. (|4]). 

Dipole branching leads to an exponential increase of 
their number with rapidity, and consequently, to an un- 
limited rise of T if formula ([6]) is applied: This would 
be inconsistent with unitarity, which requires the bound 
T < 1. It is believed that this rise is tamed by nonlinear 
effects [lij that limit the number of dipoles in the tar- 
get. Unfortunately, the latter effects have not yet been 
formulated in the framework of the dipole model. It is 
not even clear that dipoles should still be the relevant 
degrees of freedom in that regime. 

On the other hand, one knows the equation that gives 
the average scattering amplitude at rapidity Y+dY given 
the amplitude at rapidity Y, with the exact nonlinearity 
that preserves the unitarity of T. It reads: 

(T(Y + dY, r 01 )\T(Y, r 01 )) = T(Y, r 01 ) + 
d 2 r 2 rl x 



adY 



An '02' 12 



(T(Y,r 2)+T(Y,ri2)-T(Y,roi) 
-T(Y,r 02 )T(Y,r 12 )). (7) 



This equation is equivalent to the first equation in the cel- 
ebrated Balitsky hierarchy Eq. ([7]) may be obtained 
from Eq. ([2|) by averaging it over the noise v between 
rapidities Y and Y + dY. A Fourier transformation to 
position space completes the identification. In practice, 



Eq. ([7]) is useful in a regime in which T may be approxi- 
mated by its average, i.e. when (T(Y+dY, r i)\T(Y, roi)) 
may be replaced by T(Y + dY, roi), turning ([7|) into a 
closed equation for T = (T), Eq. ([3]). This is realized 
when the underlying effective number of dipoles is large, 
i.e. in the region of high density. The corresponding 
equation is the BK equation |4|, l5( and is in the universal- 
ity class of the (deterministic) FKPP equation [l5| . The 
evolution of higher order correlators of Ts that would be 
needed to be able to evolve T also outside the dense re- 
gion has still not been derived in a systematic way for 
the generic case of dipole-dipole scattering. 

Although the correct formalism to describe the evo- 
lution of T accurately in all regimes is still not avail- 
able, exact asymptotic results [16| universal enough to 
also apply to QCD amplitudes have been found from the 
approximate formulation Unfortunately, their valid- 
ity is quite reduced since the limit h\(l/a 2 ) 3> 1 had to 
be assumed, which is of course out of experimental reach. 
Hence it would be useful to understand what can be ex- 
pected quantitatively from what is known up to now, 
beyond the far asymptotics. 

The goal of this paper is to show that quite accurate re- 
sults for the evolution of the QCD scattering amplitudes 
may be extracted from a numerical study that appro- 
priately matches a Monte-Carlo simulation of the color 
dipole model, valid in the dilute regime, to the evolution 
of the amplitude given in Eq. (J7J), useful in the dense 
regime. We shall propose a factorization procedure of T, 
on an event-by-event basis, that we justify by a calcula- 
tion in the framework of a zero-dimensional model and 
motivate and test on a one-dimensional toy model. We 
then explain how this factorization could be applied to 
QCD. 



II. FACTORIZATION IN A 
ZERO-DIMENSIONAL MODEL 

In a first stage, we study a zero-dimensional model 
(transverse variables are not considered) in order to intro- 
duce the factorization rigorously. Zero-dimensional mod- 
els were investigated for possible applications to QCD 
some time ago [l7], [l8[ . Solutions for models of that kind 
have recently been discussed in Refs. QjJUll, from differ- 
ent perspectives. The solutions which have been obtained 
will help us to assess the validity of the factorization that 
we shall propose here. 



A. Definition of the model 

We investigate the time evolution of a specific Marko- 
vian process in the universality class of the zero- 
dimensional stochastic FKPP equation. We consider a 
system of n(t) particles. Between times t and t + dt, 
each particle has a probability p + — dt to split in two 
particles. For each pair of particles, there is a probabil- 
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ity dt/N that one of them is lost, and thus each given 
particle has probability p_ = (n t — l)dt/N to disappear. 
These rules completely define the process. 

There are several ways to represent this evolution. We 
may write the distribution of the number n t +dt of par- 
ticles at time t + dt given the number of particles n t at 
time t: 



P(n t+ dt\n t ) = S nt+dunt I 1 - n t dt - n f T 1 dt 



N 

+ 5n t+dt ,n t +in t dt + 6 nt+dt , nt - 1 nt ( nt N ^ dt. (8) 

This equation may be cast in the form of a stochastic 
equation for n t by first computing the mean and variance 
of nt+dt given n t , with the help of Eq. ([5]). This enables 
one to write the time evolution of nt in terms of a drift 
and of a noise of zero-mean and normalized variance, 
namely: 



dt 



N 



dn t n t (n t - 1) / n t {n t - 1) 

n t — r h V «H — vt+dt-, (9) 



N 



where v is such that (v t ) = and (vtVt 1 ) = 5(t — t 1 ). 
Note that the distribution of v depends on n t and is 
not a Gaussian. This last point is easy to understand: 
The evolution of v t is intrinsically discontinuous, since 
it stems from a rescaling of n tl which is an integer at all 
times. A Brownian evolution (i.e. with a Gaussian noise) 
would necessarily be continuous. For completeness, we 
write the statistics of vt+dt , which is easy to derive from 
the evolution of n: 



— — proba n t dt 



^t+dt 



proba 1 — n t dt 



-TF3t-% P roba 



n t (n t — l) 
N 



n t (n t -l) 
N 



dt 



dt, 



where A = n t — 



Tt t («t-1) 



and a 



N " — V n * ' N 

see well the jumps induced by the terms proportional to 
1/dt. 

Formulating the process with the help of a stochastic 
equation such as © has no real advantage at this point. 
Here, we just aimed at showing explicitely the connection 
with stochastic partial differential equations. 



i t (n t -l) 



(10) 
We 



B. Poissonian states and "Pomerons" 

Instead of looking at a state of definite occupancy, one 
could also consider the evolution of a Poissonian state 
whose occupation numbers are distributed as 



(ii) 



and follow the evolution of z t , that is to say, compute the 
probability distribution of z t +dt given z t , P{zt+dt\zt)- As 



can be easily checked, the moments of Zt+dt are the fac- 
torial moments of n t +dt- This statement may be written 
as 



E 



n-t+dt 



(n t+d t - k)l 



P(nt+dt\nt)P Zt (nt) 



dz t+ dtZt+dt p ( z t+dt\zt). (12) 



Replacing Eqs. dU) and ([TT]) in Eq. (12]), one finds 



dz t+ dtz t+ dtP( z t+dt\zt) — z t + dt 

,2 



Zt 



N 



2dt z t 



ZL 
N 



k(k - l)z. 



k-2 



(13) 



from which it is obvious that P(zt+ d t\zt) is a Gaussian 

2 2 

of mean z t + dt(z t — -^-) and variance 2dt{zt — -^). In 



the same way as one translates Eq. © into Eq. (|9]), this 
may be expressed in the form of a stochastic evolution 
equation for z t [2l| 



dz t 
dt 



zt 




] Vt+dt- 



(14) 



This is an Ito equation: r\ is a Gaussian white noise sat- 
isfying (r)t) — and (rjtrjt') = 8(t — t'). Zt is now evolving 
continuously, unlike n*. Eq. (114|) is the zero-dimensional 
version of the stochastic FKPP equation. It is of course 
consistent with Eq. ([9]), since the moments of z are the 
factorial moments of n. Starting from Eq. (|14[) and trans- 
forming it into a hierarchy for the factorial moments of 
n t , one can compute the first few orders of (n t ) in a l/N 
expansion resummed for large t 19]. The result turns 
out to be an asymptotic Borel-summable series. It reads 

MM 



A^(-l) fc+1 A- fe fc!e fc *. 



(15) 



Each term may be interpreted as the result of the evalu- 
ation of a diagram with a number fc of "Pomerons" b eing 
exchanged in the f-channel. We refer the reader to [191 ] 
for a detailed calculation in such a framework. 

The series ([15)) can be rewritten in the form of a Borel 
integral, 



(n t ) = N 1 - Ne 



which eventually may be expressed in terms of special 
functions. 

Eqs. (|T5|) and (fl~6]) are a prio ri valid for e t /N <C 1- Ac- 
tually, it was found in Ref. [19| from the exact evaluation 
of subleading orders, that its effective range of validity is 
much broader, t/N <C 1. 

The classical limit is also easily identified from 
Eq. (fl~4"l) : Indeed, a classical state has a definite value 
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FIG. 1: Bunch of dotted lines: Ten different paths for n t 
obtained from the evolution of the zero- dimensional model 
starting from a single-particle initial condition. Dashed line: 
Mean-field solution, see Eqs. (JTTJ) , d2DJ ■ Full line: Numerical 
solution for (nt) obtained by averaging over a large number 
of realizations of Eq. © . 

of z. The classical path z(t) = (nt) ~ n^ F is obtained 
by solving the deterministic part of Eq. (|14|) , namely 

j MF / MF^2 

dt ~ 4 N ■ [U) 

This approximation is expected to be valid only when 
the typical number of particles in the system is large. 
It is useful to notice that the steady fixed point of the 
evolution is (n t ) — N . 

The formulation of Eq. (|14[) is helpful for analytical 
calculations of the moments of nt- However, in order 
to describe a physical system that starts evolution with 
a definite number of particles, one needs to introduce 
complex values of the Poisson parameter z. This makes 
Eq. (|14p of little interest for numerical simulations of 
such systems. Furthermore, the analytical method for 
the computation of moments is very awkward to trans- 
pose numerically, since it leads to an asymptotic series 
that eventually needs to be resummed. Finally, it is not 
straightforward to generalize to models with spatial di- 
mensions. 



C. Solving the stochastic evolution for n: 
dense-dilute factorization 

For the difficulties mentioned above with the Poisso- 
nian state approach, we wish to take a different point of 
view on the evolution of the stochastic model. Instead 
of performing a large- N perturbative calculation directly 
of (n t ) using field-theoretical methods, we try to charac- 
terize the shape of each realization of an evolution that 
starts with no = 1 particles. What we mean by "realiza- 
tion" is a given path for n generated by the stochastic 
evolution ©. 



It is useful to visualize a few such realizations: This 
is shown in Fig. Q] for N = 5000, together with the so- 
lution to the mean field equation (|17p. One sees that 
the curves that represent n t look like the solution to the 
mean- field equation (|17[) . but with the origin of times 
translated by some random to. (The curves look also 
slightly noisy around the average trend, but the noise 
would be still much weaker for larger values of N). This 
suggests that once there are "enough" particles in the 
system, say n ^> 1, the evolution becomes essentially de- 
terministic. Hence stochasticity only manifests itself in 
the initial stages of the evolution until nt = n, but in a 
crucial way: Indeed, as seen in Fig. [TJ after averaging, 
(n t ) is significantly different from n^ F , and this differ- 
ence stems from rare realizations in which the particle 
number stays low for a long time. Therefore, in indi- 
vidual realizations, stochasticity should be taken into ac- 
count exactly as long as n t < n. Fortunately, when the 
number of particles in the system is small compared to 
the maximum number of particles N, the stochastic evo- 
lution is essentially governed by a linear equation which 
is not difficult to handle analytically. 

This heuristical discussion suggests that we may fac- 
torize the evolution in a linear stochastic evolution up to 
the time at which the number n of particles in the system 
reaches n, and continue through a nonlinear but deter- 
ministic equation, which is obtained from a mean field 
approximation to the evolution equation. As we will see, 
this simple observation leads to an elegant computation 
of (nt) that consistently agrees with the lowest order (see 
Eq. (PSD) derived in Ref. [l2|,[I!- 

Let us denote by p fi (t) the distribution of the times at 
which the number of particles in the system reaches n 
and (nt\ni) the conditional average number of particles 
at time t given that there are ni particles in the system 
at time i. One may write 

/>oo 

(nt) = / di P n(t)(n t \nt). (18) 
Jo 

So far, this expression is exact. 

We now assume that the evolution is linear when 
n t < n and deterministic for n t > n. Thus, in the pre- 
vious equation, we replace p«(t) by the solution p^ n (i) 
of the linear problem obtained by setting p + = dt and 
p_ = 0. Furthermore, we approximate (nt\ni) by the 
solution to the nonlinear evolution in the mean field ap- 
proximation (|17|) over a time interval t — i starting with 
n particles at time i, that we denote by ■ In these 

notations, we find the factorization 

POO 

(n t ) = / dtp l nt)n^ m . (19) 
Jo 

Note that for i > t (i.e. when there are less than n 
particles in the system at the considered time t), n^^- 
is like a backward evolution towards lower number of 
particles. This is not a problem since we then go back to 
the dilute regime, and the solution for (n) in that regime 
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is just obtained by taking the dilute limit of n , given 
by the solution of the equation obtained from Eq. (fT7|) 
by dropping the nonlinear term. 

We claim that this factorization is valid whenever n 
is large enough compared to 1 to justify the mean field 
approximation for the subsequent evolution, but, at the 
same time, n is small compared to N in such a way that 
the evolution up to n be linear. 

Let us express give explicit expressions for the differ- 
ent quantities that appear in Eq. P^|) . The solution to 
Eq. (JT7J) reads 



MF 



N 



1 + K e -ft-i) 



(20) 



for n -C N. It is also not difficult to show that p y ™ solves 
the equation 



•Hi) = (n - 1) / dF pl^e-^- 1 ^'). (21) 



P„ 



Its solution may be found by standard generating func- 
tion methods and takes the simple form 



p^(t) = (n-l)e- t (l- e -*) s -' 



(22) 



In the limit of interest, that is for large n and t, it boils 
down to a Gumbcl distribution 



Pt(t) = 

Replacing Eqs. flU and ([23]) in Eq. flU}, we find 



N 



dt 



ne 



1 + K e -{t-i) ' 



(23) 



(24) 



Since the Gumbel distribution is strongly damped for i < 
(by a factor e _n ), we may safely replace the lower limit 
of the integral by — oo. Finally, the change of variable 
b = -^e' - * brings the integral in the form (|16p . The result 
is independent of n, as it should be, but in numerical 
applications, one should keep in mind that the limits 1 <C 
n <C N have been assumed, and thus finite n corrections 
must be expected. 

In this picture, the Borel integral (fl6|) has a trans- 
parent interpretation: the parameter b is related to the 
factorization time i, the exponential is the distribution 
of i and the denominator corresponds to the mean-field 
part. Our method allowed to directly arrive at the lead- 
ing order result found in Refs. (l7l [l9| . without having 
to resum an asymptotic (divergent) series. 



III. A ONE-DIMENSIONAL TOY MODEL 

We shall now formulate and test this dense-dilute fac- 
torization on a one-dimensional toy model. Again, we 
consider a system made of typically N particles in its 
steady state. 



With one space dimension labeled by the real variable 
x, for large enough times, the particles form a front that 
travels towards say larger values of x when time flows 
Its position X t may be taken as the position of 
the (iV/2)-th rightmost particle. This front connects a 
dilute region, to the right, to a denser region to the left. 
There is a point in the front around which the typical 
number of objects is n. In the spirit of the previous sec- 
tion, we will treat the evolution to the right of this point 
as stochastic and solve a mean field equation to the left. 
We will explain the factorization and illustrate it numer- 
ically on a specific particle model that was introduced in 
Ref. 

When one takes a new step in time t — > t + At, the par- 
ticle at position Xi(t) is replaced, with probability At, by 
two particles at positions, where S\ and 62 are distributed 
according to ip(di)d8i (resp. ■0(^2)^2)- This rule defines 
the evolution of the number of particles in a way analo- 
gous to the dipole splitting rule in QCD. To implement a 
simple form of saturation of the number of particles, only 
the N particles which have largest positions are kept for 
subsequent evolution. 

We define T(t, x) as the fraction of these N particles 
that have positions larger than x at time t, that is 



T(t,x) 



1 N 

i=i 



(25) 



Obviously, T(t, —00) = 1 and T(t, +00) = for t large 
enough for the total number of particles in the system to 
have reached N. 

Considering for a while the N — 00 limit, the mean 
evolution of T in one elementary step in time reads 



(T(t+At,x)\T(t,x)) 

mil) I 1, (1 - At)T(t, x) + 2At J d8ip(6)T(t,x - 5) 



(26) 



This equation is the analogous of the first equation in the 
Balitsky hierarchy in the QCD case, see Eq. ([7]). In the 
infinite N limit, it is known that the large time solutions 
are traveling waves whose average velocity is determined 
by the linearized part of Eq. (|26| . If one defines 



(27) 



v{j) = In ^1 - At + 2Ai J dS e 7 V(<5)^) 



then 1^(7) admits a minimum at 7 = 70 and the large-time 
front velocity is v(jo). At finite N, the position of the 
front becomes stochastic, with nontrivial (non-Gaussian) 
statistics. The moments of the position of the front are 
known for large ./V [3] . The first two of them read 



V x t =(X t ) = v( lQ )t 
D x t =(X?) - (X t ) 2 = 



r 2 l 2 v"(lo) 
2 In 2 N 
., _ 7r 4 7ot/'(7o) 



t, 



3 In TV 



(28a) 
(28b) 
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The average position of the front was first obtained by 
considering a deterministic evolution equation with a cut- 
off [13, HH that simulates the discreteness of the particles 
in the system (24|. In our case, one would write 

T cutoS {t + At, x) = min ^1, (1 - At)T cutoS (t, x) 

+ 2At J d8^(S)T cutoS (t,x-5)j 

x 6[T cutoff (t + At, a) -1/N], (29) 

The first correction to this approximation (not shown in 
Eq. (|28ap ) is also known [la ]. It is due to particles that 
are randomly sent ahead of the deterministic part of the 
front (to the left of the cutoff), at some distance to the 
right of the tip of the front. Their multiplication through 
time evolution pulls the front forward, and generates a 
positive correction to the velocity and the dispersion in 
the front position D given in Eq. (|28b[) . One also knows 
that the system is completely renewed every t ~ 1/Z? 
units of time |23j . This remark will help the analysis of 
the numerical data. 

In our numerical implementation, we choose ip to be 
the uniform distribution in the interval [0, 1], i.e. 

i/,(g) = 6(5)6(1 - S), (30) 

and At = 0.1. Solving i>'(7o) = where v(j) is given by 
Eq. (j2"7) with these settings, we get 

70 = 1.46256- •• , 
u(7o) = 2.07006 ••• , (31) 
„"( 7o ) = 0.753472- •• . 

7o is the logarithmic slope of the falloff of the front. v(-fo) 
is the velocity in the N = oo limit. v"(7o) is a parameter 
that appears when one considers finite- N corrections. 

We first solve the complete stochastic model. We take 
an initial condition of the form X\ = • • • = Xn = (that 
is T(t = 0, x) = 6(— x)) and evolve it in time using 
the exact evolution rules. We evolve the system over 
p x 500 units in t (p is 500 in our simulation), recording 
the position of the front X ti every 8t = t% — tj_i = 500 
units of time. We compute V = (8iX t ) /St, where &iX t — 
X u -X u _ t , and D = {{8iX%) - (5iX t ) 2 )/St. The average 
is over the p periods of St units of time in one single 
realization, but for an ergodic system, averaging over the 
time evolution of one realization is like averaging over 
many independent realizations of the evolution. 

Since the system renews itself every 1/D units of time, 
and since 1/D is at most of order 100, the differences 
of successive positions SiX t are essentially independent 
random variables. Hence we can estimate the statistical 
uncertainties by splitting our set of positions 5iX t in say 
10 subsets, and by evaluating the dispersion of V and D 
between these different subsets. This helps us to provide 
error bars. The results are shown in Tab. Q] and in Figs.[H 
and |3] for N = 10 2 , 10 3 , 10 4 , 10 5 . 
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FIG. 2: Average velocity of the front V versus the equilibrium 
number of particles N for the toy model. The full curves rep- 
resent theoretical estimates. The black squares is the result 
of the mean-field calculation with a cutoff. The solution of 
the full stochastic model is shown with circles linked by large 
dashed lines. The solution of the mixed method for n = 10 
is displayed with crosses linked by short dashed lines. White 
squares denote the case n = 100, but they are almost indis- 
tinguishable from the solution of the full stochastic model. 



We also solve the deterministic evolution with a cut- 
off, defined in Eq. (|2"5|) . To this aim, we need to take 
a discretization in the x variable: We define a lattice 
with 1000 points per unit of x. The integration over S 
in Eq. (|2"9"]) is performed using the rectangle method. Al- 
though this method converges very slowly (one expects 
a systematic error of the order of 0.1% in our settings), 
it is quite well suited here because as we impose a cutoff 
at each step of the evolution, T cutoS has sharp disconti- 
nuities. The result is shown in Fig. [^together with the 
theoretical formula (|28ap . with which there is a very good 
agreement except for the lowest values of N. But that is 
expected since Eq. (|28al) was obtained in a large- N limit 
of the deterministic evolution with a cutoff. The slight 
discrepancy (of order 0.1%) visible at large N is consis- 
tently explained by the fact that our discretization and 
our integration method amount to solving a model on a 
lattice rather than the original model in the continuum, 
for which one can compute the asymptotic front velocity 
w(to), which is indeed slightly lower. 

Next, we perform a mixed evolution, applying our 
dense-dilute factorization procedure. This goes as fol- 
lows. At all times, the forward part of the front where 
T(t,x) < n/N is represented by the positions of the n 
foremost particles, and the latter are evolved using the 
exact rule (we must drop the nonlinearity, but for this 
model it only amounts to selecting the N most forward 
particles: If n is small enough, the forward particles alone 
cannot have more than N offsprings in one iteration). 
This is illustrated in Fig.[4]for a particular realization and 
for a particular choice of the parameters N and n. For 
T(t,x) > n/N, we directly evolve T using Eq. (f2"6"|) with 
the approximation (T(t + At, x)\T(t, x)) = T(t + At,x), 
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FIG. 3: Diffusion coefficient of the front D versus N. Same 
legend as in Fig. [2] 
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FIG. 4: One realization of the front for N = 100 evolved using 
the mixed method with a factorization scale n = 10. Below 
x m , the front is represented by T evolved deterministically. 
For x > x m , the position of each particle is recorded and 
evolved using the exact stochastic rules. 



forgetting about the exact position of each individual 
particle. More technically, an evolution over the time 
interval At goes as follows. We first generate particle po- 
sitions from T in the region of x between x fi — 1 and Xn- 
Indeed, in our model, it is precisely the particles present 
in that interval that may split into the forward region 
where we decided to keep the full stochasticity. This is 
related to our particular choice of ip, see Eq. (l30|) , but in 
practice, it is enough that the splittings be local enough 
in x. This is indeed the case for the fixed impact parame- 
ter BFKL evolution in QCD (see Eq. © with x = Ink 2 ). 
To get back the particle positions from the profile of T is 
extremely straightforward in this particular model: the 
position Xi of particle number i is the rightmost point for 
which T(t,Xi) — i/N holds (In more subtle models such 
as QCD, one would have to invert a relation like Eq. ([6])). 
The positions of the obtained particles are shown with 
crosses in Fig. [5] We then evolve these additional parti- 
cles together with the n forward particles using the full 
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FIG. 5: One realization of the evolution of the front of Fig. U 
(reproduced here in dashes) over the time interval At. The 
backward part is evolved through Eq. ([26} (full line). The 
crosses mark the positions of the particles that are used as 
an input for the stochastic evolution of the forward tail (they 
include the particles represented in Fig. [3] whose positions are 
tracked exactly, and additional particles generated from T). 
The circles denote the particles that are kept at time t + At: 
most of them were already present at t, and two of them have 
been produced stochastically. The forward part of T at time 
t + At is represented in short dashes. 



stochastic rule. Simultaneously, we evolve T determinis- 
tically for all x, but we cut the result at x m defined as 
the position in the front for which T(t + At, x m ) = n/N . 
Indeed, only for T > n/N do we trust the mean-field 
approximation (T(t + At, x)\T(t, x)) = T(t + At,x). In 
turn, we only keep the forward particles that have posi- 
tions larger than x m . Finally, the matching of the for- 
ward and backward parts of the front is done by requir- 
ing that T be a decreasing function of x. In the case in 
which there is a number of particles larger than n pro- 
duced ahead of x m , then T computed from the forward 
particles is larger than n/N at x m : We choose to con- 
tinue T(t + At, x m ) for x < x m until the point at which 
T(t + At,x) = T(t + At,x m ). Note that with this pre- 
scription for the matching, the particle distribution in 
the forward part of the front is exact at all times, but 
the number of particles in the backward part is slightly 
overestimated on the average. Other prescriptions may 
have been chosen: They would lead to similar quantita- 
tive results. One particular realization of the resulting 
front is shown in Fig. 

We turn to the discussion of the numerical results for 
the mixed method. Recall that the factorization scale n 
has to be at the same time larger than 1 in such a way 
that mean field evolution starting from h can be justified, 
and much smaller than N so that nonlinear effects may 
be neglected. We choose n = 10 and n — 100. This 
means that the minimum value of N for which we may 
solve the model using our method is of the order of 100 
for n = 10 and 1000 for n = 100. We follow exactly the 
same procedure as in the fully stochastic case: We evolve 
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N stochastic n = 10 n — 100 



V 


10 2 


1.8937(3) 


1.9173(5) 


— 




10 3 


1.9824(2) 


1.9910(2) 


1.9838(4) 




10 4 


2.0176(2) 


2.0202(3) 


2.0173(3) 




10 5 2.0348(2) 


2.0345(2) 


2.0331(2) 


D 


10 2 


121.2 ±3.9 


85.47 ±4.51 




xlO 3 


10 3 


41.87 ± 2.46 


33.89 ± 2.96 


44.19 ± 2.74 




10 4 


21.99 ±0.68 


16.07 ± 1.21 


18.73 ± 1.36 




10 5 


11.82 ±0.93 


8.243 ± 0.498 


10.38 ±0.68 



TABLE I: Numerical results for the one-dimensional toy 
model. The velocity V and the diffusion coefficient D of the 
front are computed in the full stochastic model and in the two 
mixed models with factorization scales n = 10 and n = 100 
respectively, for four different values of N. 



the same initial condition over p x 500 steps of time, with 
p = 500. 

The calculations of the velocity of the front are pre- 
sented in Tab. |TJ In Fig. [2j they are compared to the 
results obtained by solving the exact model. We observe 
a perfect agreement for n = 100, at the level of 0.1%. 
For n = 10, the agreement is less good: the velocity is 
overestimated by about 1% (but this means a 10% mis- 
match in the difference v (70) — V with the infinite- N 
velocity). Of course, a better agreement than that could 
not really be expected. However, part of the mismatch 
may also be related to the fact that our matching pre- 
scription between the dense and dilute regions leads to a 
slight overestimate of the number of particles in the lower 
part of the dense region, which has indeed the effect of 
increasing the front velocity. It may well be that for dif- 
ferent models which do not have the requirement that T 
be a decreasing function of x, like QCD, the agreement 
would even be better. 

The diffusion coefficient D of the front is given in Tab. [I] 
and displayed in Fig. [3] We also plot the theoretical for- 
mula (|28b[) . to guide the eye. The numerical data seem to 
agree well with the asymptotic theoretical estimate (I28bl) , 
but this must be accidental since for the considered values 
of N, corrections are expected to be huge. The diffusion 
coefficient is a quantity that is much more difficult to 
measure numerically than the velocity, since it requires 
more statistics. Our error bars are still of the order of 5 
to 10%. We see however that the full stochastic model is 
very well reproduced, within errors, by the mixed model 
with n = 100. For h = 10 instead, we notice that the 
fluctuations seem to be globally underestimated by up 
to 30% (see Tab. J]). But this should be expected since 
our method consists in replacing part of the stochastic 
evolution by a mean field approximation. 



IV. OUTLOOK FOR QCD 

We have worked on specific toy models of noisy evolu- 
tion of the FKPP type. We have proposed that, on an 
event-by-event basis, their evolution may be factorized 
into a linear stochastic part and a nonlinear determin- 
istic part. We have shown on a zero-dimensional model 
that writing down this factorization leads to a straight- 
forward computation of the average number of particles 
at leading order in t/N, without having to go through the 
resummation of a divergent series like in more standard 
approaches fl7l. [lH IT^ . 

We have conjectured that such a factorization is also 
valid for one-dimensional models. In fact, this factor- 
ization had already been implicitly assumed in previous 
works in statistical physics, essentially for the purpose of 
performing numerical calculations for very large values of 
N @, [Mil. In Ref. [3, it was even used to obtain an- 
alytical results for the front position, but only exponen- 
tially large values of N were attainable. We have shown 
here that this factorization may be extended (in practice 
numerically) to lower values of N, of the order of 100 
or 1000. The only condition is the existence of a meso- 
scopic scale of particle numbers n such that 1 <C n -C N. 
We have tested the factorization in a model that we could 
also solve numerically, and have shown that it reproduces 
quite well the exact solution even for values of N as low 
as 100. 

Of course, the dense-dilute factorization would be of no 
interest for numerical calculations at low values of N in 
cases in which the complete underlying stochastic model 
were known - systems with 100 or 1000 particles may 
often be simulated quite easily. However, as we have re- 
called in the introduction, the exact formulation of high 
energy QCD as a stochastic process (if it exists) has not 
been found yet. Instead, we have a probabilistic rule for 
the evolution of the dipole number n valid in the dilute 
regime of small amplitudes (see Eq. ([5])), which has al- 
read y b een implemented numerically by Salam [l8T ] (see 
also [27]), and a mean-field equation for evolving the am- 
plitude T in the dense regime (the Balitsky-Kovchegov 
equation 0, @, see Eq. ([7])). In our toy model, the lat- 
ter is equivalent to Eq. (|2"rJ|) and the former is like the 
evolution rule for our system of particles. In addition, 
the relationship between n and T in the dilute regime is 
needed: it is provided by Eq. © in QCD, and has its 
equivalent (Eq. (|25"|) ) in the toy model of the previous 
section. Consequently, the factorization proposed here 
should be very well suited for numerical computations of 
QCD scattering amplitudes at high energy: One should 
combine Salam's dipole Monte-Carlo and a solution of the 
Balitsky-Kovchegov equation, in a way similar to what 
we have done here for the one-dimensional toy model. 
The details will be worked out in another publication 

Extrapolating the results of our toy model study, we 
can expect to get reliable results for models in which 
the typical allowed number of particles N is larger than 
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100. This corresponds to a s < 0.1 in QCD, which is not 
far from the experimentally accessible window. Hence 
the factorization procedure outlined here may be a good 
starting point for a realistic numerical investigation of 
QCD amplitudes near the unitarity limit. Anyway, it 
is probably the best one could achieve without finding 
a realization of nonlinear saturation effects in the dipole 
model. Needless to say, the toy model that we have stud- 
ied here was tuned to make our factorization procedure as 
easy as possible to handle numerically. The many com- 
plications of QCD will make the implementation of this 
factorization a challenging issue. 

At a more theoretical level, what our study suggests 
is that universal features show up already for quite low 
values of N. This leaves us with the hope that the asymp- 
totic analytical calculations of Ref . [l6| could be extended 
to a wider range in N. 

Finally, we wish to comment on a recent proposal 
on how to address numerically the problem of QCD 
evolution at very high energies beyond the mean-field 
BK limit. In Ref. [29|, the idea that fluctuations could 
be obtained by solving a classical equation and then 
averaging over an ensemble of initial conditions was 
suggested, and was implemented numerically subse- 
quently in Ref. [3(3| ■ We note that the philosophy of this 
proposal is orthogonal to the one developed here: In 
our view, fluctuations in the saturation scale are due to 
intrinsic noise related to the discreteness of the number 
of partons, as is natural in reaction-diffusion processes. 
We do not think that the approach of Refs. [29|, yQ| 



would work for reaction-diffusion at very large times 
(i.e. asymptotic energies): For example, a universal 
distinctive feature of such processes is that the dispersion 
of the front positions (i.e. saturation scales) between 
different realizations (i.e. events) scales like VY, which 
cannot be reproduced in the approach of Refs. [29|, |30( |. 
Since our numerical method relies on the conjecture 
that QCD evolution is in the same universality class as 
reaction-diffusion processes, the two approaches would 
not lead to the same results for QCD observables (see 
Ref. |13|). However, since the statement that there is 
a correspondence between QCD and reaction-diffusion 
is still a conjecture, one has to keep open to other 
possibilities. 
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